RoutingModels

History

current version 1.6 - 30th April 2026

version	date	comment
1.0	03/Oct/2016	Original code
1.1	05/May/2023	bug corrected in RungeKutta third order
1.2	28/Nov/2023	manage reservoir when level is high
1.3	12/Dec/2023	ecological flow is defined at daily scale
1.4	24/Apr/2025	added 'noleap' when computing doy in DivertFlow
1.5	08/May/2025	bugfix writing Qout in reservoirs
1.6	30/Apr/2026	variable dt in reservoirs

License

license: GNU GPL http://www.gnu.org/licenses/

Module Description

This module includes models for discharge routing

Uses

- Units
- Utilities
- Reservoirs
- Diversions
- DataTypeSizes
- Chronos
- TableLib
- StringManipulation
- LogLib

Variables

Type	Visibility		Name
real(kind=float),	public	::	Qfilter
logical,	public	::	QfilterApply
logical,	public	::	dtAdjust	set when reservoir dt must be adjusted
integer(kind=short),	public	::	dtMin	minimum time step (s) for reservoir routing when variable step is on

Functions

private function Beta(y, B0, alpha, slope, n) result(cf)

returns the correction factor to be used in MCT algorithm (-)

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)	::	slope	bed slope (m/m)
real(kind=float),	intent(in)	::	n	manning riughness coefficient (s m^(-1/3))

Return Value real(kind=float)

private function Celerity(y, B0, alpha, slope, n) result(c)

returns the celerity (m/s)

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane [deg)
real(kind=float),	intent(in)	::	slope	bed slope (m/m)
real(kind=float),	intent(in)	::	n	manning riughness coefficient (s m^(-1/3))

Return Value real(kind=float)

private function DepthFromArea(A, B0, alpha) result(y)

returns water depth given wetted area (m) using Newton-Raphson.

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	A	area (m2)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)

Return Value real(kind=float)

private function Discharge(y, B0, alpha, slope, n) result(Q)

returns the normal discharge (m3/s) with Chezy equation.

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section(m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)	::	slope	bed slope (m/m)
real(kind=float),	intent(in)	::	n	manning riughness coefficient (s m^(-1/3))

Return Value real(kind=float)

private function NormalDepth(Q, B0, alpha, slope, n) result(nd)

returns the normal depth (m) using Newton-Raphson.

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	Q	discharge (m3/s)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)	::	slope	bed slope (m/m)
real(kind=float),	intent(in)	::	n	manning riughness coefficient (s m^(-1/3))

Return Value real(kind=float)

private function RungeKutta(hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, eFlow, freeFlow, freeFlowElevation, QoutDOY, hmin) result(stage)

solve mass continuity of reservoir using runge kutta method

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	hini	initial stage value (m)
real(kind=float),	intent(in)	::	Qin_begin	input discharge at the beginning of time step (m3/s)
real(kind=float),	intent(in)	::	Qin_end	input discharge at the end of time step (m3/s)
integer,	intent(in)	::	dt	time step (s)
type(Table),	intent(in)	::	geometry	geometry of reservoir
integer,	intent(in)	::	rk	solution order 3 or 4
integer,	intent(in)	::	typOut	type ofoutflow: 1=free flow 2=target level
real(kind=float),	intent(in)	::	hTarget	follow a target (observed) stage (m)
real(kind=float),	intent(in)	::	eFlow	environmental flow (m3/s)
real(kind=float),	intent(in)	::	freeFlow	free flow (m3/s)
real(kind=float),	intent(in)	::	freeFlowElevation	free flow elevation (m asl)
character(len=3),	intent(in)	::	QoutDOY	doy to compute Qout
real(kind=float),	intent(in)	::	hmin	minimum stage value (m)

Return Value real(kind=float)

private function TopWidth(y, B0, alpha) result(b)

returns the topwidth of a triangular, rectangular or trapezoidal cross section

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)

Return Value real(kind=float)

private function Velocity(y, B0, alpha, slope, n) result(v)

returns the normal velocity (m/s) with Chezy equation

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)	::	slope	bed slope (m/m)
real(kind=float),	intent(in)	::	n	manning riughness coefficient (s m^(-1/3))

Return Value real(kind=float)

private function WettedArea(y, B0, alpha) result(a)

returns the wetted area of a triangular, rectangular or trapezoidal cross section

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)

Return Value real(kind=float)

private function WettedPerimeter(y, B0, alpha) result(p)

returns the wetted perimeter of a triangular, rectangular or trapezoidal cross section

Arguments

Type	Intent		Name
real(kind=float),	intent(in)	::	y	water stage (m)
real(kind=float),	intent(in)	::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)	::	alpha	angle formed by dykes over a horizontal plane (deg)

Return Value real(kind=float)

Subroutines

public subroutine DivertFlow(time, div, Qin, Qout)

Discharge routing through diversion. Discharge is split between river and diversion channel

Arguments

Type	Intent		Name
type(DateTime),	intent(in)	::	time
type(Diversion),	intent(inout)	::	div
real(kind=float),	intent(in)	::	Qin	Input discharge at time t + dt (m3/s)
real(kind=float),	intent(out)	::	Qout	downstream discharge at time t + dt (m3/s)

public subroutine LevelPool(time, dt, dtrk, Qin_t, Qin_tplusdt, pool)

Discharge routing through reservoir. Given Qin, find Qout and stage References:

Arguments

Type	Intent		Name
type(DateTime),	intent(in)	::	time
integer,	intent(in)	::	dt	main time step (s)
integer,	intent(in)	::	dtrk	time step to run runge-kutta (s)
real(kind=float),	intent(in)	::	Qin_t	Input discharge at time t (m3/s)
real(kind=float),	intent(in)	::	Qin_tplusdt	Input discharge at time t + dt (m3/s)
type(Reservoir),	intent(inout)	::	pool

public subroutine MuskingumCunge(dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

Solve Muskingum-Cunge equation.

Arguments

Type	Intent	Optional		Name
integer(kind=short),	intent(in)		::	dt	time increment (s)
real(kind=float),	intent(in)		::	dx	reach length (m)
real(kind=float),	intent(in)		::	Qin	input discharge at current time (m3/s)
real(kind=float),	intent(in)		::	Pin	input and output at previous time (m3/s)
real(kind=float),	intent(in)		::	Pout	input and output at previous time (m3/s)
real(kind=float),	intent(in)		::	Qlat	lateral flow at current and previous time (m3/s)
real(kind=float),	intent(in)		::	Plat	lateral flow at current and previous time (m3/s)
real(kind=float),	intent(in)		::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)		::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)		::	slope	bed slope (m/m)
real(kind=float),	intent(in)		::	n	manning roughness coefficient (s m^(-1/3))
real(kind=float),	intent(out)		::	Qout	output discharge at current time
real(kind=float),	intent(out)		::	depth
real(kind=float),	intent(out)		::	width
real(kind=float),	intent(in),	optional	::	k	flood wave travel time used to compute C1, C2, and C3 (s)
real(kind=float),	intent(in),	optional	::	x	used to compute C1, C2, and C3 (-)

public subroutine MuskingumCungeTodini(dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

Solve Muskingum-Cunge-Todini equation. When parameters k and x are passed, C1, C2, and C3 are computed using the passed values. When x = 0, linear reservoir model. When x = 0.5 and k = dt translation only

Arguments

Type	Intent	Optional		Name
integer(kind=short),	intent(in)		::	dt	time increment (s)
real(kind=float),	intent(in)		::	dx	reach length (m)
real(kind=float),	intent(in)		::	Qin	input discharge at current time (m3/s)
real(kind=float),	intent(in)		::	Pin	input and output at previous time (m3/s)
real(kind=float),	intent(in)		::	Pout	input and output at previous time (m3/s)
real(kind=float),	intent(in)		::	Qlat	lateral flow at current and previous time (m3/s)
real(kind=float),	intent(in)		::	Plat	lateral flow at current and previous time (m3/s)
real(kind=float),	intent(in)		::	B0	bottom width, = 0 for triangular section (m)
real(kind=float),	intent(in)		::	alpha	angle formed by dykes over a horizontal plane (deg)
real(kind=float),	intent(in)		::	slope	bed slope (m/m)
real(kind=float),	intent(in)		::	n	manning roughness coefficient (s m^(-1/3))
real(kind=float),	intent(out)		::	Qout	output discharge at current time
real(kind=float),	intent(out)		::	depth
real(kind=float),	intent(out)		::	width
real(kind=float),	intent(in),	optional	::	k	flood wave travel time used to compute C1, C2, and C3 (s)
real(kind=float),	intent(in),	optional	::	x	used to compute C1, C2, and C3 (-)

RoutingModels Module

Contents

Variables

Functions

Subroutines

History

License

Module Description

Uses

Variables

Functions

private function Beta(y, B0, alpha, slope, n) result(cf)

Arguments

Return Value real(kind=float)

private function Celerity(y, B0, alpha, slope, n) result(c)

Arguments

Return Value real(kind=float)

private function DepthFromArea(A, B0, alpha) result(y)

Arguments

Return Value real(kind=float)

private function Discharge(y, B0, alpha, slope, n) result(Q)

Arguments

Return Value real(kind=float)

private function NormalDepth(Q, B0, alpha, slope, n) result(nd)

Arguments

Return Value real(kind=float)

private function RungeKutta(hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, eFlow, freeFlow, freeFlowElevation, QoutDOY, hmin) result(stage)

Arguments

Return Value real(kind=float)

private function TopWidth(y, B0, alpha) result(b)

Arguments

Return Value real(kind=float)

private function Velocity(y, B0, alpha, slope, n) result(v)

Arguments

Return Value real(kind=float)

private function WettedArea(y, B0, alpha) result(a)

Arguments

Return Value real(kind=float)

private function WettedPerimeter(y, B0, alpha) result(p)

Arguments

Return Value real(kind=float)

Subroutines

public subroutine DivertFlow(time, div, Qin, Qout)

Arguments

public subroutine LevelPool(time, dt, dtrk, Qin_t, Qin_tplusdt, pool)

Arguments

public subroutine MuskingumCunge(dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

Arguments

public subroutine MuskingumCungeTodini(dt, dx, Qin, Pin, Pout, Qlat, Plat, B0, alpha, slope, n, Qout, depth, width, k, x)

Arguments